
extern double interpol_2dl(const double *x, const double *y, const double *u, double xp, double yp)
     {
      double a, b, c;

      a = ((u[0]*y[1]-u[1]*y[0])*(x[1]*y[2]-x[2]*y[1]) - (u[1]*y[2]-u[2]*y[1])*(x[0]*y[1]-x[1]*y[0]))\
             /((y[1]-y[0])*(x[1]*y[2]-x[2]*y[1]) - (y[2]-y[1])*(x[0]*y[1]-x[1]*y[0]));

      b = ((u[0]-u[1])*(y[1]-y[2])-(u[1]-u[2])*(y[0]-y[1]))/((x[0]-x[1])*(y[1]-y[2])-(x[1]-x[2])*(y[0]-y[1]));

      c = ((u[0]-u[1])*(x[1]-x[2])-(u[1]-u[2])*(x[0]-x[1]))/((y[0]-y[1])*(x[1]-x[2])-(y[1]-y[2])*(x[0]-x[1]));

      return a + b*xp + c*yp;
     }
